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ABSTRACT 


In this thesis we present an algorithm for the precise determination of the Mie 
extinction efficiency parameter. The mathematical representation of the Mie parameters is 
in the form of an infinite series, and any technique that could be found to accelerate the 
convergence of the Mie series would have great commercial and military application. 
Results are presented that show the comparison of the rate of convergence obtained by 
directly summing the individual terms of the extinction efficiency parameter and the rate 
obtained using an existing series acceleration technique. It was found that the acceleration 
method we employed, known as the Levin method of series transformation, proved 
unsuccessful in accelerating the convergence of the Mie series. However, other 


acceleration techniques exist and should be explored. 
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I. INTRODUCTION 


Mie scattering is the scattering of electromagnetic radiation by primarily spherical 
particles whose diameters are comparable to the wavelength of the incident radiation 
[Ref. 1: p. 303]. In nature, we see this manifested as the white appearance of clouds. Due 
to the distribution of particle sizes in the clouds, the cumulative effect of the Mie scattered 
light waves on the water droplets is that the whole spectrum of scattered radiation 
combines to make the cloud appear white. In other scenarios which involve the direct 
use of modern technology, Mie scattering has many applications where it is of vital 
importance to have a highly efficient and reliable algorithm for signal processing and for 
analyzing the properties of Mie scattered waves. A few examples include the reflection of 
a radar signal off of a cloud of dispersed particles and the absorption and scattering of 
radiation from soot particles in the atmosphere [Ref. 2]; Mie scattering as a technique for 
the sizing of air bubbles [Ref. 3]; and for remote sensing applications. Further applications 
can be found in fields as diverse as astrophysics [Refs. 4 - 5], physical chemistry [Refs. 6- 
7], and a unique style of painting that employs the Mie scattering phenomenon {Ref. 8]. 

In order to determine the amplitude, extinction efficiency, or scattering efficiency 
of an electromagnetic wave that has undergone Mie scattering, it is necessary to compute 
the mathematical representation of these quantities, or Mie parameters, each of which is 
analytically represented by an infinite sum. For the purpose of this thesis research, our 


attention has been focused on calculating the Mie parameter known as the extinction 


efficiency factor, or Q.x:, which describes the total effect of scattering and absorption in 
removing radiation from the incident beam. 

As will be shown later in Chapter IV, the number of terms required for the infinite 
series to converge numerically to a given accuracy is directly proportional to the ratio 
between the wavelength of the incident radiation and the circumference of the scattering 
particle. The thrust of this thesis research has been to determine if this infinite series can 
be caused to converge in dramatically fewer iterations through the utilization of a series 
acceleration technique. The basic idea is to transform the Mie series into another series 
that converges faster. While numerous series acceleration methods are known, one of the 
most powerful is known as the Levin method [Ref. 9: p. 35]. The ultimate goal of this 
thesis research has been to determine if the Levin method can significantly accelerate the 
convergence of the Mie series. 

As an example of the power of the Levin method (detailed in chapter 3), consider 


the alternating harmonic series that represents the natural logarithm of 2, namely, 





es sj n+] 
in(2)= 5°! m1-p tant (1) 


The convergence of this series is extremely slow. By summing the individual terms of 
Eq. (1), over a million and a half terms must be taken to achieve an accuracy of six digits: 
1,565,239 terms to be exact. The C code used to determine this value and the program 
output are included as Appendices A and B, respectively. By employing the Levin 
method, however, the same degree of numerical accuracy can be reached after only six 


iterations! The ANSI C code for the Levin method operating on the alternating harmonic 


series is included as Appendix C. After the fifteenth term of the Levin method is reached, 
the accepted value of accuracy out to the sixteenth digit is obtained [Ref. 10: p. 113]. 

This is an astounding reduction of computer processing time and is an amazing 
demonstration of the power of the Levin method of series acceleration. (We note that, to 
obtain eight digit accuracy by straightforward summation of the series would require over 
one hundred million terms; nine digit accuracy was effectively unattainable using the 
UNIX system resources available for this research.) Table 1 shows a comparison of the 
partial sum obtained from Eq. (1) and those calculated with the Levin method. Asa 
further comparison, the plot shown in Fig. 1 was produced to show the relative speed (i.e., 
number of iterations required) with which the Levin method and the brute force method of 


summing the individual terms of the alternating harmonic series converge to 99 


In2 = 0.6931471805599453 


Partial Sum 
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0.725372 0.693 1471805599453 


Table 1. Levin Method on Logarithm of 2. 
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Figure 1. Levin Method vs. Infinite Series for In(2) 


percent of the correct value of In(2). Appendices D through F show the pertinent ANSI C 
and Matlab code. Note that the Levin method required only three terms to reach 99 
percent of the value of In(2), while summing the individual terms of the series required 
sixty-five terms. As a final comparison, the time required to reach six digits of accuracy 
using the brute force method was determined (Appendices A and B), and the code used in 
Appendix C was altered to determine the number of times the Levin method could reach 
the same degree of accuracy in the same time span (11.7 seconds). The value obtained 
was approximately 20,000, indicating that in this instance the Levin method is about 
20,000 times faster than summing the individual terms of the infinite series. 

In light of the preceding example, the question arises as to whether the Levin 
method can accelerate the Mie series. A computer program was written to calculate the 
individual terms of the infinite series representation of the Mie parameter Q-,, and to 
determine the number of these terms required for the series to converge to an accuracy of 
one part in 10°. Next, the Levin method was applied to this infinite series, and the number 
of iterations of the Levin method required to achieve a similar accuracy was likewise 
determined. A comparison of these two results showed that the Levin method did not 
succeed in accelerating the convergence of the Mie series. In fact, the number of 
iterations required for the Levin method to converge was comparable to that obtained by 
simply summing the infinite series representation of the Mie parameter. Considering the 
algorithm which the Levin method utilizes, this method when applied to the Mie series will 


actually result in increased processor time. Although the Levin method proved 


unsuccessful in accelerating the convergence of the Mie series, further research in this area 
is still needed as there are other series acceleration methods worth exploring 
[Ref. 11: p. 56]. 

In the next section, we discuss the relevant Mie parameters and the infinite series 
for computing the extinction coefficient. In Chapter III we discuss the Levin convergence 
acceleration method. In Chapter IV we discuss our results for applying the Levin method 


to the Mie series. 


Il. THE MIE SCATTERING PARAMETERS 


A. BACKGROUND 

There are four basic Mie scattering parameters. The first two pertain to the 
scattering amplitude and describe the complex amplitudes of the perpendicular and parallel 
components (with respect to the scattering plane) of the electric field vector. The third 
Mie parameter is the scattering efficiency factor, while the fourth parameter is the 
extinction efficiency factor. All four of these expressions are given in dimensionless terms, 
and the evaluation of these four basic functions relies on the proper evaluation of four 
coefficients: the Mie coefficients a, and b,, and the angular coefficients z, and 7. 

These four expressions can be described using three basic parameters: a 
dimensionless size x of the scattering sphere, an index of refraction m (usually complex), 
and a scattering angle @ (relative to the forward direction). The dimensionless size 


parameter x is given by 


wg (22, 7 circumference of sphere | (2) 
A wavelength 


where k =27/A is the wave number, r is the radius of the sphere, and / is the 

wavelength of the incident radiation. The complex index of refraction m is given by 
m=v-Ki , (3) 

where v is the real part of the index of refraction and « is the imaginary part. A complex 


refractive index indicates an absorbing sphere, such as soot. A purely real index of 


refraction is indicative of a non-absorbing sphere; for the purpose of this thesis research, 
only cases of purely real refractive index were considered. 
B. SCATTERING AMPLITUDE PARAMETERS 

Let the electric field vector amplitude of the scattered radiation field be given by 
Asc. This radiation field can in turn be expressed in terms of the scalar perpendicular and 
parallel components A, and A2. We can define dimensionless, complex amplitudes S, and 
S2 by multiplying the amplitudes A; and A, by the free-space propagation constant k, 


which in turn can be represented by an infinite converging series [Ref. 12: p. 13]: 








kA, = S\(m,x,0)= fa,z,(u)+b,7,(u)} (4a) 
— l 
kA, = S,(m,x,) = 2 . : 1) {b,7,(u)+4,7,(u)} (4b) 


The Mie coefficients a, and b, are functions of the index of refraction m and the 
size x, while the angular coefficients z, and 7, are functions of sz@cos@ only. The latter 
coefficients are defined in terms of Legendre polynomials and their derivatives. The 
coefficients a, and b, are functions of spherical Bessel functions of the first, second, and 
third kind, and can be expressed in a variety of ways. 

C. SCATTERING EFFICIENCY AND EXTINCTION PARAMETERS 
It can be shown that the differential scattering cross section for unit incident flux is 
do(m,x,0)=4 Asc: Agc(m,x,A)do , (5) 


where dq@ is an element of the solid angle. 


Following Diermendjian [Ref. 12: p. 13], by representing the unpolarized incident 
radiation as the sum of two independent and linearly polarized components of equal flux, 


the scattering cross section can be expressed as 
Gen. X i= eos x, 0) = 3h A'+A,4,)do , (6) 


where (2 = 477 1s the solid angle. 
The scattering efficiency factor Q,.(m,x) is obtained by normalizing o, by the 
sphere’s geometrical cross section 27°: 


Om Ge = Fal) “ ; 





+|b,|°) ° (7) 








n=] 
The total extinction cross section and efficiency factor, which includes the 
contribution due to absorption, can be similarly defined. This leads to the following 


expression for extinction efficiency factor: 
2 N 

Onn =z DA (2n+1)Re(a,+5,) . (8) 
n=l 


As noted earlier for the case of a purely real index of refraction m, there will be no 
absorption, and thus the extinction efficiency factor Q-x; will be identically equal to the 
scattering efficiency factor Osea . 

D. THE MIE COEFFICIENTS 
The simplest and most elegant expressions for the Mie coefficients are as follows 


[Ref. 13: p. 195]: 


ec 


a, = er iey (kr) ’ (9a) 


» - nr) — 9), (7) 
”"— krh_,(kr)—nh (kr) 


(9b) 
where jn(kr) is a spherical Bessel function of the first kind with order n, h’,(kr) is a Bessel 
function of the third kind, known as a Hankel function, and is defined as 
h nkr) = j,(kr) +iy, (Ar), (10) 
and y,(kr) is a spherical Bessel function of the second kind with order n. 
Another common form of the Mie coefficients is that adopted by van de Hulst 


[Ref. 14: p. 123]: 


_ AnlZ)V n(x) mV) (11a) 
"A, (2)G (x) — me (x) 


p, = AZ a(2)= WAC) ais 


mA,(2)6,,(x) — C3, (x) 


where z = mx , while y, and ¢, are Ricatti-Bessel functions which are defined as 


follows: 
V (x)= xi,(x) , | (12a) 
LAAE—xy,@) (12b) 
C(x) = V(x) tz, (x) - (12c) 


Here, j, and y, are spherical Bessel functions of the first and second kind, respectively. 
For computational purposes, it is best to express the Mie soc in a form 

conducive to separating their real and imaginary components. Through the use of 

recursion formulas and circular functions [Ref. 12: p. 16-19], the Mie coefficients may be 


written in the alternative form [Ref. 15: pp. 16-17]: 








(“6 az) LY ACK) — Va) 
a, = a ; (13a) 
= Deh, ()— Cn. (x) 
XxX 
jm) + y (x)-y,\(X) 
bu= x (13b) 


{m+ (2) 6.) 


These forms of the Mie coefficients were programmed, using the C language, to determine 
the Mie coefficients and hence the extinction efficiency factor Qex . 
oe CALCULATING @Q,, USING ANSI C 

The author wishes to thank the writers of "Numerical Recipes in C" [Ref. 16] for 
enabling the circumvention of the reinvention of the wheel, in that their code for 
computing ordinary and spherical Bessel functions and their derivatives was utilized in this 
thesis research for computing the Mie coefficients and calculating the extinction efficiency 
factor. Of noteworthy interest is that said code cites Jerry Lentz of the Naval 
Postgraduate School as the author of an improved technique for calculating Bessel 
functions through the use of continued fractions. 

The code used for calculating the ordinary Bessel functions is incorporated into a 
structure which returns the values of the Bessel functions of the first and second kind, and 
their first derivatives. The spherical Bessel function code simply generates a normalization 
factor, makes a call to the function which calculates ordinary Bessel functions (of half-odd 
integer order), then in like fashion returns values of the spherical Bessel functions of the 


first and second kind, and their first derivatives. The code for calculating the ordinary 


Bessel functions also makes calls to two additional functions, "beschb" and "chebev", also 
obtained from "Numerical Recipes in C" [Ref. 16]. The “beschb” function is used to 
evaluate the gamma functions present in the Bessel function formula, calling the “chebev” 
function in the process to perform Chebyshev expansion. 

Once the interfacing codes for calculating ordinary and spherical Bessel functions 
were brought together, a program was written to assemble the components of the Mie 
coefficients as expressed in Eqs. (13a) and (13b). Again, the Mie coefficients are 
functions of the dimensionless size parameter x and the index of refraction m. 

Next, a loop was created that calculated the individual terms of the extinction 
efficiency factor (one of the Mie parameters), each term of which relied on a calculated 
value of the Mie coefficients. As is the case of all of the Mie parameters, the extinction 
efficiency factor 1s represented as an infinite sum. Thus, the loop for calculating the 
efficiency extinction factor was summed until the value converged to one part in 10°. 
With each increment of the loop counter, the Bessel functions associated with the Mie 
coefficients were of subsequent higher order, 1.e., for n=1, J,(x) would be a Bessel 
function of the first order; for n=2, J,(x) would be a Bessel function of the second order, 
and so on. 

As a test for the accuracy of the code and the correctness of the programming, a 
routine was created that calculated the extinction efficiency factor for size parameter x 
varying from 0 to 9, using a step size of 0.01 and applying 6 different values for the 
refractive index. A plot showing all six curves was produced from the data obtained from 


this program and is shown in Fig. 2. This plot very closely resembles that which was 


12 


produced by Van de Hulst [Ref. 14: p. 151], with the exception that Van de Hulst's plot 
was intentionally smoother as a result of not showing fine detail. The C code used to 
produce data for Fig. 2 is included as Appendix G, with the output data being sent to a file 


called “Qout”. The Matlab code used to produce Figure 2 is included as Appendix H. 





ii. THE LEVIN METHOD 


A. SERIES ACCELERATION METHODS 

The main idea behind the Levin method, as in any series acceleration technique, is 
to transform a given slowly converging series into another series that converges faster. 
We note here the distinction between a sequence and a series. An infinite sequence (S,) is 
an unending progression of numbers S, which may be real or complex. A sequence 
converges if a number, S, exists so that, corresponding to every positive number €, no 


matter how small, a number n, can be found such that 





S,-S|<e for n>n,. In this case, 
the sequence (S,,) is said to converge to the limit S as ” tends to infinity. An infinite series, 
on the other hand, is the sum of an infinite sequence. Let u,,u,,...,u,,... be an infinite 
sequence of numbers, real or complex. Let the sum uw, +u,+---:u, be denoted by S,, 
which is called the n” partial sum. Then, if the sequence of partial sums (S,) converges to 
a limit, S, the infinite series u,+u,+--- 1s said to be convergent, or to converge to the 
sum S. The connection with the Levin method is to transform the sequence of partial 
sums into another sequence that converges faster, 1.e., requires fewer terms to converge to 
the limit. 

Given a sequence of real or complex numbers (S,,) which converges to S, this 
sequence can be transformed into another sequence (7;,). A trivial example of such a 


sequence transformation is 


i= UN. (14) 


To be useful, however, the following properties for the transformed sequence must hold 
[Ref. 11: p. 1]: 
1. (7,) must converge, 


2. (7,) must converge to the same limit as (S,,), and 


3. (T,) must converge to S faster than (S,,); that is, lim(T, - S)/(S, -—S)=0. 


If property 3 holds, the transformation T is said to accelerate the convergence of the 
sequence (S,,), or that the sequence (7,,) converges faster than (S,). An example of a 


sequence transformation meeting these three conditions is Aitken’s A’ process [Ref. 11: 





pp. 1-7]: 
AS AS 
TS See). So 15 
n As. n ( Se n+} ( ) 
or, 
2 
P =D Pat Osea (16) 
See _ 28 1 a; Ss 


where A is the difference operator defined by Av, =v,,,—v, and A‘*'y, = A‘y_,, —A‘v.. 


The A’ in the denominator of Eq. (15) accounts for the process’s name. 

The Levin transformation, as with most other transformation algorithms, is a 
particular case of what is known as the E-transformation, which is the most general 
sequence transformation. Transformations belonging to the E-transformation class include 
[Ref. 11: p. 56]: Richardson polynomial extrapolation, Shanks’ transformation, the Levin 


transformation, and many other historically known acceleration methods. We note that 


series acceleration is an active field of research in numerical analysis, with the Levin 
method and E-transformation having been discovered only in the past 25 years. 
B. THE E-ALGORITHM 

Following Brezinski and Zaglia [Ref. 11: pp. 56-57], the E-transformation is based 
on the following relation: 

DO ae) ae (17) 
where S,, is an element of the sequence to be transformed, the a’s are unknown scalars, the 
gi(n)’s are given auxiliary sequences (which may depend on terms of the sequence S,, 
itself), and where k is a fixed integer. Rewriting Eq. (17) in terms of S,, 

S,=S+a,g,(n)+-:-+a,2,(n) . (18) 
The basic idea behind Eq. (18) is to attempt to fit the actual behavior of S,, as a function 
of n, so that it may be extrapolated smoothly to the (unknown) limit, S. This fitting is 
achieved through the particular choices of the functions g,(m). By incrementing 7 in Eq. 
(18) to n + k, one has k + 1 equations in k + 1 unknowns (the limit S and the scalars 


ae 7” ,a,). By solving these equations, we obtain a sequence of estimates for the limit, 


which we denote by E’”: 


SS. oe Sa 
Zi) <-> Sita k) 
Bp =) elt (19) 
gin): g(n+k) 
g,(n) +: g(ntk) 


It is assumed that the determinant in the denominator of Eq. (19) is not equal to zero. 
Finally, note that Eq. (18) is the kernel of the transformation Eq. (19); that is, Eq. (18) is 


the set of sequences for which there exists the sequence S such that, for all n, E - =. 


The E-algorithm is the recursive algorithm that allows one to compute the numbers 


E{” without actually computing the determinants in Eq. (19). Given the following rules, 


ja ee | (20a) 


i te) nH=@)1,... and i=1,2,... , (20b) 
then for k= 1,2,...andn=0,1,..., the main rule for the E-algorithm is 


) (n) Dy Ale (”) 
(nm) _ n k-1 elle Gg 
Ey 7 Ee (ntl) _ An) Sk-1,k ? (21) 


k-1,k k-1,k 


where the gi”), ’s are auxiliary sequences computed by the following auxiliary rule: 


Gai. = eae 

(n) _ (rn) k-l,i kK-l,t(n) an 

Eki = Aa = (n+) (n) *Sk-Lk 5 i= k ate l, k +2, eco (22) 
Si-1k ~ Sk-1,k 


C. LEVIN’S TRANSFORMS 

Levin’s method of generating non-linear transformations for increasing the rate of 
convergence of sequences [Ref. 17: pp. 371-388] was first introduced in 1973. For 
Levin’s generalized transformation, the auxiliary sequences denoted by the gi({)’s in the 


E-transformation given in Eq. (19) take on the particular form: 


x TAS 
2,1) == eee (5) 


n 


where x, and y, are themselves auxiliary sequences. 


Levin’s transforms are basically generalizations of Aitken’s A’ process and of the 
E-transformation corresponding to the first column of the E-algorithm as shown in 
Eg. (21) [Ref. 11: p. 113], denoted as E,. The kernel of Aitken’s process is 
S,-S=a-AS, , (24) 
while the kernel of the transformation E, is of the form 
S,-S=t,-g(n) . (25) 


By definition in the Levin method, g(7) is taken to be an arbitrary polynomial of degree 


(k—1) of the quantity (n). Also, t, denotes the n” element of the series (as for example 


the n’” term of the Mie series). Thus, the sequences of Eq. (25) can be expressed as 


S,-S=t, a, +d, (n) +-+++a, (ny (26) 


Following Brezinski and Zaglia, multiplying both sides of Eq. (26) by (n) and 


rearranging terms yields 


(n\” “a = a,(n)*” +a, (n) +e--+a, . (27) 


By applying the operator A‘ to both sides of Eq. (27), the right hand side becomes 


identically zero, since it is a polynomial of degree (k-—1) inn. Therefore, for all n, 


ato 5.5) =(i (28) 


n 


Since A‘ is a linear operator, the above expression can be also be expressed as 


a(n &.) =S-A eo) (29) 


ni "i 
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The numbers E‘”), denoted now as the Levin estimate L\”) , are therefore given by 


(n) _ a"((n)" ‘S, / a 


(30) 
AM) 
Using the well-known formula [Ref. 11: p. 115] 
N - k sk 
ANU, = (-l) Cyt (31) 
k=0 
where C;, is the binomial coefficient defined by 
N 
= oct Malian : (32) 
x} kN Bb)! 
Eq. (30) then becomes 
= k n-15 
(1) (n+ ky Sak /e(N—H) 
ie = es ; (33) 
S'(-1)'(n+ bw /k (N—k) 


= 
tl 
— 


n+k 
Equation (33) defines the sequence of Levin estimates for the limit of the series. 
D. ENCODING THE LEVIN TRANSFORM 
In the instances where the Levin code (Appendices C and L) was utilized, Ty was 
used to represent the partial sum S_,, from Eq. (33) above. The sum 7y is related to the 
sequence tf; by the following expression: 
a 
=e (34) 
= 


] 


In other words, ¢, represents the individual term for the loop which calculates the 


extinction efficiency factor Q,,, and 7,, is the value for which Q,_ converges to one part 
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in 10°. Equation (33) was split into two individually calculated terms, P[N | and Q| N | : 
representing the numerator and denominator, respectively. Consequently, the following 


equations result: 


k= wp 


— 


PUN ]= (1) (e)"" vt few —4) , and ( 35) 


Q[N]= 5 (-1)'(k)"" as jk (N-k)! . (36) 


k=1 ly 


Thus, Eq. (30) becomes 
i= ea (37) 


Equation (34) was encoded in order to calculate the Levin estimate of Q,., and the 


number of terms needed for the Levin transform to cause the infinite series representation 


of Q,,, to converge was subsequently determined. This number was compared with that 


which was obtained by direct summing of the individual terms of the infinite series 


representation of Q,,, to convergence in order to determine if the Levin method 


accelerated the convergence. 
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Figure 2. Extinction Efficiency Factor vs. Size Parameter 
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IV. RESULTS 


After completing the Mie code and producing values for the extinction efficiency 
factor that agreed with the tabulated results of Van de Hulst, the result obtained from 
applying the Levin method of series acceleration to the infinite series representation of a 
Mie scattering parameter was somewhat anti-climactic, albeit not entirely unanticipated. 
Simply put, the Levin method did not accelerate the convergence of the extinction 
efficiency factor. In point of fact, the number of iterations required for the Levin method 
to converge to a given accuracy was approximately the same, less one or two terms, as the 
number of terms required by direct summation. Size parameters ranging from 0 to 100 
were used in making this determination. Since the Levin method requires calls to 
additional functions that calculate the log of a gamma function, the log ofa factorial, and a 
binomial coefficient at each iteration of the summation loop — the codes for which were 
obtained from “Numerical Recipes in C” [Ref. 16] — this method is actually slower than 
direct summation. It should also be noted that, for size parameters x > 108, the Levin 
method fails entirely because the individual terms f, become vanishingly small, producing 
undefined output. The code used to determine Q,., via summing the infinite series and the 
code utilizing the Levin estimate are included as Appendices I and J, respectively. 

The C code used to produce the data used for the plot shown in Figure 2 was 
slightly altered to determine the number of iterations required for the Mie series to 
converge as a function of size parameter. The result is shown in Figure 3. Since both 


methods of convergence required the same number of iterations, the code used for 
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summing the individual terms of the Mie series (Appendix G) was chosen to generate the 
data used for the plot shown in Figure 3. A linear relationship can easily be seen for all 
indices of refraction, with the number of terms in the series proportional to the size 
parameter. 

While a tremendous amount of effort was put forth to come to the conclusion that 
the Levin method is unsuccessful in accelerating the convergence of the Mie series, the 
work was not all for naught. An important question was answered, leading the way for 


further research in this area. 


a =N EN aN 
‘o] NO +e 0) 


# iterations to converge 
© 





Size parameter x 


Figure 3. Rate Convergence vs. Size Parameter 


24 


V. CONCLUSION 


In this thesis research we assembled the code to calculate the Mie coefficients, 
integral components of all four of the Mie scattering parameters. Next, we calculated the 
Mie extinction efficiency factor by summing the terms of the infinite series representation 
to convergence. Several different indices of refraction were used, and a plot of the 
extinction efficiency factor versus size parameter was produced. These results were found 
to agree closely with that of Van de Hulst’s authoritative work, thereby ensuring the 
accuracy of the code. Finally, an algorithm for the Levin method of series transformation 
was incorporated into the existing code. The rate of convergence was subsequently 
determined and compared with the rate of convergence achieved by summing the 
individual terms of the Mie series. Based on the results, we arrived at several important 
conclusions. 

First, the number of terms required to converge was found to be linearly related to 
the size parameter. This included size parameters both within and outside the Mie regime. 
The linear relationship continued until a limitation was reached, which leads us to our 
second conclusion. 

The magnitude of the size parameter was found to impose a limit on the extent to 
which the Levin method was able to cause the Mie series to converge. While a size 
parameter of 108 is well beyond the Mie scattering region and within the optical region, it 


still points to a computational limit of the Levin method. 


2s 


Most importantly, we have concluded that the Levin method of series 
transformation did not accelerate the convergence of the Mie series. Several simulations 
were run using size parameters ranging from less than one up to one hundred, and indices 
of refraction ranging from one to two. All results showed a rate of convergence, insofar 
as the number of iterations required, equal to the rate obtained by summing the individual 
terms of the Mie series. This is an important result. There are other series transformation 
methods in existence that may accelerate the convergence of the Mie series, and these 


should be explored in their own right. 
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APPENDIX A: ALTERNATING HARMONIC SERIES CODE IN ANSI C 


/* Thesis Program to calculate number of terms for alternating harmonic series to reach progressive */ 
/* digits of accuracy */ 

/* LT Brian Johnson */ 

/* Compiler: Borland C++ Ver. 5.0 */ 

/* File Name: |n2.c */ 


# include <stdio.h> 
# include <math.h> 
# include <time.h> 


int n; 
int places = 1; 


double AltHarmSum = 0.0; 
double RoundedNew = 0.0; 
double RoundedOld = 0.0; 


float accuracy = 10.0; 


main() 


{ 


n= 0; 


do { 
n+ |; 
AltHarmSum+=(double)( pow(-1,n+1)/n); 
RoundedNew = floor(AltHarmSum *accuracy+.5)/accuracy; 


if ( RoundedNew == RoundedOld ) 
{ 
printf("\nIt took %od terms to reach %d digits of accuracy:", n, places); 
printf("\n Sum = %f", AltHarmSum); 
printf("\n Time elapsed: %3.2f sec.\n", clock()/le6 ); 
places+=1; 
accuracy*=10.0; 


} 
RoundedOld = RoundedNew; 


} while (places <= 15); 
} 


Ze 
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APPENDIX B. OUTPUT OF ALTERNATING HARMONIC SERIES CODE 


It took 12 terms to reach 1 digits of accuracy: 
Sum = 0.653211 
Time elapsed: 0.00 sec. 


It took 271 terms to reach 2 digits of accuracy: 
Sum = 0.694989 
Time elapsed: 0.00 sec. 


It took 1417 terms to reach 3 digits of accuracy: 
Sum = 0.693500 
Time elapsed: 0.01 sec. 


It took 177341 terms to reach 4 digits of accuracy: 
Sum = 0.693150 
Time elapsed: 1.33 sec. 


It took 229300 terms to reach 5 digits of accuracy: 
Sum = 0.693145 
Time elapsed: 1.71 sec. 


It took 1565239 terms to reach 6 digits of accuracy: 


Sum = 0.693147 
Time elapsed: 11.70 sec. 
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APPENDIX C: LEVIN ALGORITHM ON ALTERNATING HARMONIC SERIES 


/* Thesis Program to apply Levin method to alternating harmonic series and compute number of terms */ 
/* for series to converge */ 

/* LT Brian Johnson */ 

/* Compiler: Borland C++ Ver. 5.0 */ 

/* File Name: In2lev.c */ 


#include <stdio.h> 

#include <stdlib.h> 

#include <stddef.h> 

#include <math.h> 

#define MAX 15 

int k; 

double T[MAX+1], ttMAX+1]; 
double a; 


main() 


T[0] = t[0] = 0.0; 
k=1; 


printf("\n\ni T[i] ti] P[i}/Q[i}\n"); 
while(k<=MAX) 


{ 
t{k] = (double)( pow(-1,k+1)/k); 
T[k] = T[k-1] + t[k]; 
k=k+1; 


a = levin(); 


} 


double levin() 
{ 


double PIMAX+1],Q[MAX+1]; 
int 1; 
l=1; 
P[0] = Q[0] = 0.0; 
while(I<=MAX) 
P{1]=P{I-1]+pow(-1,1)*pow(,MAX-1)*T[I]A[I] *bico(MAX,]); 
Q{l]}]=Q[1-1]+pow(-1,1)*pow(,MAX-1)/t[]]*bico(MAX,]); 
l=1+1; 
printf("\n%d %f “of %f", i, TLi], tli}, Plij/Q[i]}); 


} 
_ (P[MAX]/Q[MAX]); 
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double bico(int n, int k) 
{ double factin(int n); 
return floor(0.5+exp(factln(n)-factin(k)-factin(n-k))); 


; 
double factln(int n) 
{ double gammin(float xx); 
static float a[ 101]; 
if (n <= 1) return 0.0; 
if (n <= 100) return a[{n] ? af{n] : (a{fn]=gammln(n+1.0)); 
else return gammIn(n+1.0); 
3 
double gammIn(float xx) 
{ double x,y,tmp,ser; 


static double cof[6]={76. 18009172947 146, -86.50532032941677, 
24.01409824083091,-1.231739572450155, 
0.1208650973866179e-2,-0.5395239384953e-5}; 

int j; 

Y~XXX, 

tmp=x+5.5; 

tmp = (x+0.5)*log(tmp); 

ser=1.000000000190015; 

for (j=0;j<=5;j++) ser += coffj]/++y; 

return -tmptlog(2.50662827463 10005 *ser/x); 

3 
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APPENDIX D: ANSI C CODE TO PRODUCE DATA FOR ALTERNATING 
HARMONIC SERIES PLOT 


/* Thesis Program to produce data for plot of convergence of alternating harmonic series to In(2) */ 
/* LT Brian Johnson */ 

/* Compiler: Borland C++ Ver. 5.0 */ 

/* File Name: In2plot.c */ 


# include <stdio.h> 
# include <math.h> 
# include <time.h> 


int n; 
int places = 1; 


double AltHarmSum = 0.0; 
double RoundedNew = 0.0; 
double RoundedOld = 0.0; 


float accuracy = 10.0; 


main() 
{ 
oman seline<= 51; nt) { 
AltHarmSum+=(double)( pow(-1,n+1)/n); 
RoundedNew = floor(AltHarmSum*accuracy+.5)/accuracy; 


printf("\n%d %1.3f", n, log(2) - fabs(log(2) - AltHarmSum)); 


} 
} 
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APPENDIX E. ANSI C CODE TO PRODUCE DATA FOR PLOT OF LEVIN 
METHOD ON ALTERNATING HARMONIC SERIES 


/* Thesis Program to produce data for plot of Levin method on the convergence of alternating */ 
/* harmonic series to In(2) */ 

/* LT Brian Johnson */ 

/* Compiler: Borland C++ Ver. 5.0 */ 

/* File Name: In2levplot.c */ 


#include <stdio.h> 
#include <stdlib.h> 
#include <stddef.h> 
#include <math.h> 
#include "nrutil.h" 
#include "nrutil.c" 
#include <time.h> 


#define MAX 16 


double levin(); 

double bico(int, int); 

double factln(int); 

double gammIn(float); 

double T[MAX+1], t(IMAX+1]; 
double a; 


int k, 1; 
main() 


T[0] = t[0] = 0.0; 
k=1; 
while(k<=MAX) 


{ 
t[k] = (double)( pow(-1,k+1)/k); 
T[k] = T[k-1] + t[k]; 
Keatcta: 

} /* end while */ 


a = levin(); 
\ /* end main */ 
double levin() 
{ 
double PIMAX+1], QUMAX+1]; 
double commonTerm; 
int N; 
RO Cl] = 0.0; 


for (N = 1; N <= MAX; N++) 
35 


{ 
for (k = 1; k <= N; k++) 
{ 
commonTerm = pow(-1,k)*pow(k,N-1)/A[k]*bico(N,k); 
P{k] = P{k-1] + commonTerm* T{k]; 
Qk] = Q[k-1] + commonTerm; 
} 


printi(\n%d %32.311f", k-1, log(2) - fabs(log(2) - P[N]/Q[N])); 
} /* end for N loop */ 


printf("\n\n"): 
return (P[MAX]/Q[MAX]); 
} 


double bico(int n, int k) 
{ double factln(int n); 
return floor(0.5+exp(factln(n)-factln(k)-factln(n-k))); 


3 
double factln(int n) 
{ double gammIn(float xx); 
Static float a[101]; 
if (n <= 1) return 0.0; 
if (n <= 100) return a[n] ? a[n] : (a[n]J=gammin(n+1.0)); 
else return gammin(n+1.0); 
} 
double gammin(float xx) 
{ double x,y,tmp,ser; 


static double cof[6]={76. 18009172947146,-86.50532032941677 
24.01409824083091,-1.23 1739572450155, 
0.1208650973 866 179e-2,-0.5395239384953e-5}; 

int J; 

Y~X>XX, 

tmp=x+5.5; 

tmp = (x+0.5)*log(tmp); 

ser=1.000000000190015; 

for §=0;j<=5;j++) ser += cof]j]/+-4y; 

return -tmptlog(2.50662827463 10005*ser/x); 

3 


? 
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APPENDIX F. MATLAB CODE TO PRODUCE PLOT OF CONVERGENCE OF 
ALTERNATING HARMONIC SERIES 


/* Thesis Program to produce data for plot of Levin method on the convergence of alternating */ 
/* harmonic series to In(2) */ 

/* LT Brian Johnson */ 

/* Software: Matlab Ver. 4.2C */ 

/* File Name: In2comp.m */ 


load In2graf -ascil 
load In2levgr -ascii 
x = In2graf(:,1); 

y = In2graf(:,2); 


x2 = In2lever(:,1); 
y2 = In2levegr(:,2); 


B = ones(1,71); 

figure(1) 

C = 0.99*B; 

plot(x-1,y/log(2),x-1,C,':', x2-1,y2/log(2)) 
xlabel('Number of iterations’) 

ylabel(‘Partial Sum / 1n2') 

title (‘Levin Method vs Infinite Sum for In(2)') 
axis({0 70 0 1]) 

gtext('0.99 -->') 
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APPENDIX G. ANSI C CODE TO PRODUCE DATA FOR EXTINCTION 
EFFICIENCY FACTOR AND CONVERGENCE RATE VS SIZE PARAMETER 


/* Thesis Program to produce data for plotting OQ, vs size parameter */ 


/* LT Brian Johnson */ 
/* Compiler: Borland C++ Ver. 5.0 */ 
/* File Name: Miegraph.c */ 


#include <stdio.h> 
#include <stdlib.h> 
#include <stddef.h> 

/* #include <iostream.h> */ 
#include <math.h> 
#include “nrutil.h" 
#include "nrutil.c" 

#define NRANSI 

#define EPS 1.0e-10 
#define FPMIN 1.0e-30 
#define MAXIT 10000 
#define XMIN 2.0 

#define PI 3.141592653589793 
#define RTPIO2 1.2533141 
#define NUSE1 5 

#define NUSE2 5 


struct bessVec { 
double rj, ry, rjp, ryp; 
He 


struct complex { 
double real, imag; 
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void assignBessels(struct bess Vec* besselsPtr, double r)dummy, double rydummy, double rjpdummy, 
double rypdummy) 


besselsPtr -> rj = rjdummy; 
besselsPtr -> ry = rydummy; 
besselsPtr -> rjp = rjpdummy; 
besselsPtr -> ryp = rypdummy; 


} 

struct bessVec sphbes(double, int); /* Function Prototype */ 

struct bessVec bessjy(double, double); /* Function Prototype */ 

Struct bessVec sphBessels_n; /* j's, y's & deriv's for Ricatti-Bes */ 
struct bessVec bessels_nPlusHalf; /* for calculating J_n +/- half */ 
struct bessVec bessels_ nMinusHalf; 

struct complex a_n, b_n; /* Mie coefficients */ 

struct complex zeta_n; /* Ricatti-Bessel fn */ 

struct complex zeta_nMinus1; 

struct complex denom; /* denominator in Mie coefficients */ 
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double x; 

double m; 

double J nMinusHalf, J nPlusHalf: 
double A_ factor; 

double psi_nMinus1]; 

double psi_n; 

double chi_n; 

double nOverx; 

double aBraces, bBraces; 

double numer; 

double realDenom; 

double sum_QextOld, sum_Qext; 
double Q_ext[900][6]; 

int num[900][6]; 

double conv = 1.0e-8; 


int 1; 
int J; 


int n; 


main() 


/* size parameter */ 

/* refractive index (non-absorbing spheres) */ 
/* for calculating A_factor */ 

/* psi-prime divided by psi */ 

/* Ricatti-Bessel fn 1st kind, initialized */ 


/* Ricatti-Bessel fn 2nd kind */ 


/* numerator in Mie coefficients */ 

/* rationalized denominator */ 

/* Extinction efficiency factor */ 

/* array of Q ext's, fn's of m and x */ 

/* number of iterations to converge, f(x,m) */ 


/* counter for x loop */ 
/* counter for mArray loop */ 
/* # iterations to converge, also order of bessel function */ 


1 
double mArray[6] = {1.25, 1.33, 1.44, 1.50, 1.55, 2.0}; 


for = 0; j <= 5; j++) { 


m = mArrtay[j]; 


for(x = 0.01; x <=9.0; x += 0.01) { 


n=0; 


/* initialize n for each x */ 


bessels_nPlusHalf = bessjy(m*x, 0.5); /* initialize for A_factor */ 
J_nPlusHalf = bessels_nPlusHalf.1; 


psi_n = sin(x); 
zeta_n.real = psi_n; 
zeta_n.tmag = cos(x); 


sum Qext = 0.0; 


do { 
n=n+l; 
J_nMinusHalf = J_nPlusHalf: 
psi_nMinus1 = psi_n; 
zeta_nMinus]1 = zeta_n; 
sum_QextOld = sum_Qext; 


/* initializing zeta */ 


bessels_nPlusHalf = bessjy(m*x, n + 0.5); /* calculating A factor */ 
J_nPlusHalf = bessels_nPlusHalf.n; 
A_factor = (J_nMinusHalf/J_nPlusHalf) - (n/(m*x)); 
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sphBessels_n = sphbes(x, n); /* j's and y's for Ricatti Bess fn's */ 
psi_n = x * sphBessels_n.yj; 

chi_n= -x * sphBessels_n.ry; 

zeta_n.real = psi_n; 

zeta_n.imag = chi_n; 


nOverx = n/x; 
aBraces = (A_factor /m + nOverx); 
bBraces = (A_factor * m + nOverx); 


numer = aBraces*psi_n - psi_nMinus1; 
denom.real = (aBraces * zeta_n.real) - zeta_nMinus1.real; 
denom.imag = (aBraces * zeta_n.imag) - zeta_nMinus1.imag; 


realDenom = (denom.real*denom. real + denom.imag*denom.imag); 


a_n.real = (numer * denom.real)/realDenom; 
a_n.imag = (numer * denom.imag)/realDenom; 


numer = bBraces*psi_n - psi_nMinus1; 
denom.real = (oBraces * zeta_n.real) - zeta_nMinus1.real; 
denom.imag = (bBraces * zeta_n.imag) - zeta_nMinus].imag; 


realDenom = (denom.real*denom.real + denom.imag*denom. imag); 


b_n.real = (numer * denom.real)/realDenom; 
b_n.imag = (numer * denom.imag)/realDenom; 


sum_Qext += (2*n + 1)*(a_n.real + b_n.real); 
} while (fabs(sum_Qext - sum_QextOld) > conv); /* end do-while loop */ 


Q ext[(int)(x*100-1)][j] = 2.0/(x*x)*sum_Qext; 
num|[(int)(x*100-1)][j] = n; 


; /* end inner x-loop */ 
} /* end outer j-loop (for various m values) */ 


for (1 = 0; 1 < 899; i++) { 
printf("\n%f ", (i+1.0)/100.0); 


for G = 0; j <= 5; j++) 
printf(of ", Q extfilf); 
/* — printf("%d_", num{i}{j]); */ 


} /* end i-j loop */ 


return 0; 
; 


struct bessVec sphbes(double x, int n) 


if 
« 


4] 


struct bessVec tempsphbessels; 
struct bess Vec bessels; 


double sj, sy, sjp, Syp; 
void nrerror(char error_text[]); 
double factor, order; 


if (n < 0 || x <= 0.0) nrerror("bad arguments in sphbes"); 
order=n+0.5; 


bessels = bessjy(x,order); 
factor=RTPIO2/sqrt(x); 
sj=factor*bessels.1; 
sy=factor*bessels.ry; 
sjp=factor*bessels.rjp-(sj)/(2.0*x); 
syp=factor*bessels.ryp-(sy)/(2.0*x); 


assignBessels(&tempsphbessels, sj, sy, sjp, SYP); 


return tempsphbessels; 
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struct bess Vec bessjy(double x, double xnu) 
{ 


void nrerror(char error_text[]); 


struct bess Vec tempBessels; 


double rj, ry, rp, ryp; 


void beschb(double x, double *gam1, double *gam2, double *gampl, 
double *gammi); 

int i,isign,} nl; 

double a,b, br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli, dr, ¢,f,fact,fact2, 
fact3 ,ff, gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,r1, 
111 ,rjmu,npl ,rjpl,rjtemp,ry 1 ,rymu,rymup,rytemp,sum,sum1, 
temp,w,x2,x1,x12,xmu,xmu2; 


if (x <= 0.0 || xnu < 0.0) nrerror("bad arguments in bessjy"); 
nl=(x < XMIN ? (int)(xnut+0.5) : IMAX(0,(int)(xnu-x+1.5))); 
xmu=xnu-nl; 
xmu2=xmu*¥xmu; 
x1=1.0/x; 
xi2=2.0* x1; 
w=x12/PI; 
isign=1; 
h=xnu* xi; 
if (h < FPMIN) h=FPMIN; 
b=x12*¥ xnu; 
d=0.0; 

=h; 
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for (=1;1<=MAXIT;i++) { 


} 


b+= x12; 

d=b-d; 

if (fabs(d) < FPMIN) d=FPMIN; 
c=b-1.0/c; 

if (fabs(c) < FPMIN) c=FPMIN; 
d=1.0/d; 

del=c*d; 

h=del*h; 

if (d < 0.0) isign = -isign; 

if (fabs(del-1.0) < EPS) break; 


if (i > MAXIT) nrerror("x too large in bessjy; try asymptotic expansion"); 
rl=isign*FPMIN; 
npl=h* Hl; 


nll=rl; 


rjpl=rpl; 
fact=xnu* x1; 
for (=nl;I>=1;1--) { 


rtemp=fact*rl+1jpl; 
fact = x1; 

rpl=fact* rjtemp-rl; 
rl=rjtemp; 


$ 

if (jl == 0.0) rjl=EPS: 
feypl/rl; 

if (x < XMIN) { 


x2=0.5*x; 
pimu=PI*xmu; 
fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimu)); 
d = -log(x2); 
e=xmu*d; 
fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e); 
beschb(xmu,&gam1,&gam2,&gampl,&gammi); 
ff=2.0/PI*fact*(gam1*cosh(e)+gam2* fact2*d); 
c=exp(e); 
p=e/(gampl*PI); 
q=1.0/(e*PI* gammi); 
pimu2=0.5*pimu; 
fact3 = (fabs(pimu2) < EPS ? 1.0 : sin(pimu2)/pimu2); 
r=PI* pimu2* fact3 *fact3; 
c=1.0; 
d = -x2*x?2; 
sum=ff+1*q; 
sum1=p; 
for G=1;i<=MAXIT;1++) { 

ff=(i* ff+ p+q)/G*i-xmuz2); 

c *= (d/l); 

p /= (i-xmu); 

q /= (i+xmu); 

del=c*(ff+1r*q); 

sum += del; 

dell=c* p-i*del; 
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} else { 


sum! += dell; 
if (fabs(del) < (1.0+fabs(sum))*EPS) break: 
j 
if (i > MAXITT) nrerror("bessy series failed to converge"); 
rymu = -sum, 
ryl = -sum1*x12; 
rymup=xmu*xi*rymu-ry1!; 
rymu=w/(rymup-f*rymu); 


a=0.25-xmu?2; 
p = -0.5*x1: 
q=1.0; 
br=2.0*x; 
bi=2.0; 
fact=a*xi/(p*p+q*q); 
cr=br+q* fact; 
ci=bi+p* fact; 
den=br*br+bi*bi; 
dir=br/den; 
di = -bi/den; 
dlr=cr* dr-ci* di; 
dli=cr*di+ci*dr, 
temp=p*dlr-q*dli; 
g=p*dlit+q* dlr; 
p=temp; | 
for (=2;i<=MAXIT;1++) { 
a += 2¥*(i-1); 
bi += 2.0; 
dr=a* dr+br,; 
di=a*di+bi; 
if (fabs(dr)+fabs(di) < FPMIN) dr=FPMIN; 
fact=a/(cr*cr+ci*ci); 
cr=br+cr* fact; 
cl=bi-ci* fact; 
if (fabs(cr)+fabs(ci) < FPMIN) cr=FPMIN; 
den=dr* dr+di* di; 
dr /= den; 
di /= -den; 
dlr=cr*dr-ci* di; 
dli=cr* di+ci*dr, 
temp=p*dlr-q*dli; 
g=p*dlit+q*dlr; 
p-temp, 
if (fabs(dir-1.0)+fabs(dli) < EPS) break; 
} 
if (i > MAXIT) nrerror("cf2 failed in bessjy"); 
gam=(p-f)/q; 
rymu=sqrt(w/((p-f)*gam+q)); 
ymu=SIGN(qmu,q)); 
rymu=rjmu* gam; 
rymup=rymu*(p+q/gam); 
ry l=xmu*xi*rymu-rymup; 
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fact=rjmu/rjl; 

rj=rjll *fact; 

rp=np! *fact; 

for G=1;1<=nl;i++) { 
rytemp=(xmuti)*xi2* ry l-rymu; 
rymu=ry 1; 
ry 1=rytemp; 

; 

ry-rymu, 

ryp=xnu*xi*rymu-ry1; 


assignBessels(&tempBessels, rj, ry, Np, ryp); 


return tempBessels; 


void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi) 


{ 


} 


double chebev(double a, double b, double c{], int m, double x); 

double xx; 

static double cl{] = { 
-1.142022680371172e0,6.516511267076e-3, 
3.087090 17308e-4,-3.470626964e-6,6.943764e-9, 
3.6780e-1 1,-1.36e-13}; 

static double c2{] = { 
1.843740587300906e0,-0.076852840844786e0, 
1.271927 136655e-3,-4.971736704e-6,-3.3126120e-8, 
2.423 10e-10,-1.70e-13,-1.0e-15}: 


Xx=8.0*x*x-1.0; 

* gam 1=chebev(-1.0,1.0,c1, NUSE1,xx); 
*gam2-chebev(- 1.0, 1.0,c2, NUSE2,xx); 
*sampl= *gam2-x*(*gam1); 

*oammi= *gam2+x*(*gam1); 


double chebev(double a, double b, double c[], int m, double x) 


{ 
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void nrerror(char error_text[]}); 
double d = 0.0, dd = 0.0, sv, y, y2; 
int j; 


if ((x-a)*(x-b) > 0.0) nrerror("x not in range in routine chebev"); 
y2 = 2.0*(y=(2.0*x-a-b)/(b-a)); 
for j=m-1;j>= 1; j-) { 

sv = d; 

d = y2*d - dd + c[j]; 

dd = sv; 


} 
return y*d - dd + 0.5*c[0]: 


#undef EPS 
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#undef FPMIN 
#undef MAXIT 
#undef XMIN 
#undef PI 
#undef NRANSI 
#undef NUSE1 
#undef NUSE2 
#undef RTPIO2 
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APPENDIX H. MATLAB CODE TO PLOT OUTPUT OF CODE LISTED 
IN APPENDIX G. 


load Qout -ascil 
x = Qout(:,1); 
yl = Qout(:,2); 
y2 = Qout(:,3); 
y3 = Qout(:,4); 
y4 = Qout(:,5); 
y5 = Qoutt(:,6); 
y6 = Qout(:,7); 


step = 0.25; 

xinc = step:step:8.75; 

for i = 1:(length(xinc)-1) 
ylinc(i) = y1(100*xinc(i)+1); 
y2ince(i) = y2(100*xinc(i)+1); 
y3inc(i) = y3(100*xinc(i)+1); 
y4inc(i) = y4(100*xinc(i)+1); 
ySince(1) = y5(100*xinc(i)+1); 
y6inc(1) = y6(100*xinc(i)+1); 

end 


plot(x,y1,x,y2,x,y3,x,v4,x,y5,x,y6, xinc,ylinc,'.’, xinc,y2inc,'x', xinc, y3inc,'*', xinc,y4inc,'+’, 
xinc,ySinc,'o’, xinc, y6inc,'.') 

axis({0,9,0,6]) 

xlabel('Size parameter x'), ylabel(‘Qext') 

gtext('m = 2') 

gtext('m = 1.55') 

gtext(‘m = 1.50') 

gtext('m = 1.44’) 

gtext('m = 1.33') 

gtext('m = 1.25') 


load Qoutn -ascli 
x = Qoutn(:,1); 
yl = Qoutn(:,2); 
y2 = Qoutn(:,3); 
y3 = Qoutn(:,4); 
y4 = Qoutn(:,5); 
y5 = Qoutn(:,6); 
y6 = Qoutn(:,7); 


step = 0.25; 
xinc = step:step:8.75; 


for i= 1:(length(xinc)-1) 
y linc(i) = y1(100*xinc(i)+1); 
y2inc(i) = y2(100*xinc(i)+1); 
y3inc(i) = y3(100*xinc(1)+1); 
y4inc(1) = y4(100*xinc(i)+1); 
y5inc(i) = y5(100* xinc(i)+1); 


47 


yoinc(i) = y6(100* xinc(i)+1); 
end 


plot(x,y1,'k',x,y2,'m',x,y3,'c'.x,y4,'1'.x,y5,'g',x,y6,'b’, xinc,y linc,'.k’, xinc,y2inc,'xm', xinc, y3inc,'*c’, 
xinc,y4inc,'tr', xinc,y5inc,'og', xinc, y6inc,'.b’) 


xlabel('Size parameter x'), ylabel('# iterations to converge’) 
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APPENDIX I. ANSI C CODE FOR CALCULATING EXTINCTION 
EFFICIENCY FACTOR 


/* Thesis Program for calculating Q.,, as a function of user-defined size parameter */ 


/* LT Brian Johnson */ 
/* Compiler: Borland C++ Ver. 5.0 */ 
/* File Name: mie.c */ 


#include <stdio.h> 
#include <stdlib.h> 
#include <stddef.h> 
#include <math.h> 
#include “nrutil.h" 
#include "nrutil.c" 

#define NRANSI 

#define EPS 1.0e-10 
#define FPMIN 1.0e-30 
#define MAXIT 10000 
#define XMIN 2.0 

#define PI 3.141592653589793 
#define RTPIO2 1.2533141 
#define NUSE1 5 

#define NUSE2 5 


struct bessVec { 


double rj, ry, rjp, ryp; 
2 


struct complex { 
double real, imag; 
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void assignBessels(struct bessVec* besselsPtr, double rj}dummy, double rydummy, double rjpdummy, 
double rypdummy) 


besselsPtr -> rj = rjdummy; 
besselsPtr -> ry = rydummy; 
besselsPtr -> rjp = rjpdummy; 
besselsPtr -> ryp = rypdummy; 


} 
struct bessVec sphbes(double, int); /* Function Prototype */ 
struct bessVec bessjy(double, double); /* Function Prototype */ 


struct bessVec sphBessels_n; 

struct bessVec bessels_nPlusHalf; 

struct bessVec bessels_nMinusHalf; 

struct complex a_n, b_n; /* Mie coefficients */ 

struct complex zeta_n; 

struct complex zeta _nMinus]; 

struct complex denom; /* denominator in Mie coefficients */ 
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double x = 109.0; /* size parameter */ 

double m = 2.0; /* refractive index (non-absorbing spheres) */ 
double J. nMinusHalf, J nPlusHalf; 

double A_factor; 


double psi_nMinus!; /* Ricatti-Bessel fn Ist kind, initialized */ 
double psi_n; 

double chi_n; /* Ricatti-Bessel fn 2nd kind */ 

double nOverx; 

double aBraces, bBraces; 

double numer; /* numerator in Mie coefficients */ 
double realDenom: /* rationalized denominator */ 


double sum_QextOld, sum_Qext, Q ext; /* Extinction efficiency factor */ 
double conv = 1.0e-8; 


int n= 0; /* loop counter, also order of bessel function */ 
main() 

{ 

bessels_nPlusHalf = bessjy(m*x, 0.5); /* initialize for A_factor */ 


J_nPlusHalf = bessels_nPlusHalf.1; 


psi_n = sin(x); 
zeta_n.real = psi_n; /* initializing zeta */ 
zeta_n.imag = cos(x); 


sum_Qext = 0.0; 


do { 
n=n+t]; 
J_nMinusHalf = J_nPlusHalf; 
psi_nMinus 1 = psi_n; 
zeta_nMinus] = zeta_n; 
sum_QextOld = sum_Qext; 


bessels_nPlusHalf = bessjy(m*x, n + 0.5); /* calculating A factor */ 
J_nPlusHalf = bessels_nPlusHalf.rj; 
A_factor = (J_nMinusHalf/J_nPlusHalf) - (n/(m*x)); 


sphBessels_n = sphbes(x, n); /* j's and y's for Ricatti Bess fn's */ 
psi_n = x * sphBessels_n.1j; 

chi_n= -x * sphBessels_n.ry; 

zeta_n.real = psi_n; 

zeta_n.imag = chi_n; 


nOverx = n/x; 
aBraces = (A_factor / m + nOverx); 
bBraces = (A_factor * m + nOverx); 


numer = aBraces*psi_n - psi_nMinus1; 


denom.real = (aBraces * zeta_n.real) - zeta_nMinus1.real; 
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denom.imag = (aBraces * zeta_n.imag) - zeta_nMinus]1.imag; 
realDenom = (denom.real*denom.real + denom.imag*denom.imag): 


a_n.real = (numer * denom.real)/realDenom; 
a_n.imag = (numer * denom.imag)/realDenom; 


numer = bBraces*psi_n - psi_nMinus]1; 


denom.real = (bBraces * zeta_n.real) - zeta_nMinus1.real; 
denom.imag = (bBraces * zeta_n.imag) - zeta_nMinus1.imag; 
realDenom = (denom.real*denom.real + denom.imag* denom. imag); 


b_n.real = (numer * denom.real)/realDenom; 
b_n.imag = (numer * denom.imag)/realDenom; 


sum_Qext += (2*n + 1)*(a_n.real + b_n.real); 
} while (fabs(sum_Qext - sum_QextOld) > conv); 
Q ext = 2.0/(x*x)*¥sum_Qext; 


printf("\n\nThe Extinction Efficiency Factor is: Qext = %f", Q ext, "\n\n"); 
printf("\nNumber of iterations: %d\n", n); 

printf("\nx = %3.2f m = %1.2f \n\n", x, m); 

return 0; 

} 


struct bessVec sphbes(double x, int n) 

{ 
struct bessVec tempsphbessels; 
struct bess Vec bessels; 


double sj, sy, sjp, syp; 
void nrerror(char error_text[]); 
double factor, order; 


if (n < 0 || x <= 0.0) nrerror("bad arguments in sphbes"); 
order=n+0.5; 


bessels = bessjy(x,order); 
factor=RTPIO2/sqrt(x); 
sj=factor*bessels. rj; 
sy=factor*bessels.ry; 
sjp=factor*bessels.rjp-(sj)/(2.0*x); 
syp=factor*bessels. ryp-(sy)/(2.0*x); 


assignBessels(&tempsphbessels, sj, sy, sjp, Syp); 


return tempsphbessels; 
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struct bessVec bessjy(double x, double xnu) 


t 
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void nrerror(char error_text[]); 


struct bess Vec tempBessels; 


double rj, ry, np, ryp; 


void beschb(double x, double *gam1, double *gam2, double *gampl, 
double *gammi); 

int 1,1sign,],n1; 

double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2 , 
fact3 ff, gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,1,11, 
rll ,rymu,npl ,rjpl,rjtemp,ry 1 ,rymu,rymup,rytemp,sum,sum1, 
temp,w,x2,x1,x12,xmu,xmu2; 


if (x <= 0.0 || xnu < 0.0) nrerror("bad arguments in bessjy"); 
nl=(x < XMIN ? (int)(xnu+0.5) : IMAX(0,Gnt)(xnu-x+1.5))); 
xmu=xnu-nl; 
xmu2=xmu*xmu; 
x1=1.0/x; 
x12=2.0* x1; 
w=x12/PI, 
isign=1; 
h=xnu*x1; 
if (h < FPMIN) h=FPMIN; 
b=x12* xnu; 
d=0.0; 
c=h; 
for G=1;1<=MAXIT;i++) { 
b += xi2; 
d=b-d:; 
if (fabs(d) < FPMIN) d=FPMIN; 
c=b-1.0/c; 
if (fabs(c) < FPMIN) c=FPMIN; 
d=1.0/d; 
del=c*d; 
h=del*h; 
if (d < 0.0) isign = -isign; 
if (fabs(del-1.0) < EPS) break; 
} 
if (i > MAXIT) nrerror("x too large in bessjy; try asymptotic expansion"); 
rl=isign*FPMIN; 
npl=h* Hl; 
nll=n); 
ypl=npl; 
fact=xnu*xi; 
for J=nl;1>=1,1--) { 
rtemp=fact*rl+rpl; 
fact = x1; 
rpl=fact*ntemp-rl; 
rl=rjtemp; 
} 
if (jl == 0.0) Hl=EPS; 
a2 


f=rnpl/rl; 
if (x < XMIN) { 


} else { 


x2=0.5*x; 
pimu=PI*xmu; 
fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimuy)); 
d = -log(x2); 
e=xmu*d; 
fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e); 
beschb(xmu,&gam1,&gam2,&gampl,&gammi); 
ff=2.0/PI*fact*(gam 1*cosh(e)+gam2*fact2*d); 
e=exp(e); 
p=e/(gampI* PI); 
q=1.0/(e*PI* gammi); 
pimu2=0.5*pimu; 
fact3 = (fabs(pimu2) < EPS ? 1.0 : sin(pimu2)/pimu2); 
r=PI*pimu2* fact3*fact3; 
c=1.0; 
d = -x2*¥x2; 
sum=ff+r* q; 
sum 1=p; 
for G=1;1<=MAAXIT;1++) { 

ff=(* ff+ p+q)/(i*i-xmuz2); 

c *= (d/i); 

p /= (i-xmu); 

q /= (i+xmu); 

del=c*(ff+r*q); 

sum += del; 

dell=c*p-i*del; 

sum] += dell; 

if (fabs(del) < (1.0+fabs(sum))*EPS) break; 
$ 
if (i > MAXIT) nrerror(“bessy series failed to converge”); 
rymu = -sum, 
ry1 = -sum1*xi2; 
rymup=xmu*xi*rymu-ry 1; 
rymu=w/(rymup-f*rymu); 


a=0.25-xmu2; 

p = -0.5*x1; 
q=1.0; 

br=2.0*x; 

bi=2.0; 
fact=a*xi/(p*p+q*q); 
cr=br+q* fact; 
ci=bitp* fact; 
den=br*br+bi* bi; 
dr=br/den; 

di = -bi/den; 
dir=cr*dr-ci* di; 
dli=cr*dit+ci*dr,; 
temp=p*dlr-q*dli; 
q=p*dlitq*dlr; 
p=temp; 
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for G=2;1<=MAXIT;1+-+) { 
a += 2*(1-1); 
bi += 2.0: 
dr=a*dr+br; 
di=a*dit+bi; 
if (fabs(dr)+fabs(di) < FPMIN) dr=FPMIN; 
fact=a/(cr*cr+ci*c1); 
cr=brtcr* fact; 
ci=bi-ci* fact; 
if (fabs(cr)+fabs(ci) < FPMIN) cr=FPMIN; 
den=dr*dr+di*di; 
dr /= den; 
di /= -den; 
dir=cr*dr-ci* di; 
dli=cr*di+ci*dr; 
temp=p*dlr-q* dh; 
q=p*dlitq*dlr, 
p=temp; 
if (fabs(dlr-1.0)+fabs(dli) < EPS) break; 
} 
if (@ > MAXIT) nrerror("cf2 failed in bessjy"); 
gam=(p-f)/q; 
rymu=sqrt(w/((p-f)*gamt+q)); 
rymu=SIGN(ymu,y)); 
rymu=rjmu* gam; 
rymup=rymu*(ptq/gam); 
ry 1=xmu*xi*rymu-rymup; 
} 
fact=ryymu/r1; 
rj=ry11*fact; 
np=np! *fact, 
for (Q=1;1<=nl;i++) { 
rytemp=(xmuti)*x12* ry 1 -rymu; 
rymu=ry1; 
ry1=rytemp; 
3 
ry=rymu, 
ryp=xnu*xi*rymu-ry1; 


assignBessels(&tempBessels, rj, ry, mp, Typ); 


return tempBessels; 


void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi) 


{ 


double chebev(double a, double b, double c[], int m, double x); 

double xx; 

static double cl[] = { 
-1.142022680371172e0,6.516511267076e-3, 
3.08709017308¢e-4,-3.470626964e-6,6.943 764e-9, 
3.6780e-11,-1.36e-13}; 

static double c2[] = { 
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1.843740587300906c0,-0.076852840844786e0, 
1.271927136655e-3,-4.971736704e-6,-3.3 126 120e-8, 
2.423 10e-10,-1.70e-13,-1.0e-15}; 


XX=8.0*x*x-1.0; 
*gam1=chebev(-1.0,1.0,c1, NUSE1,xx); 
* gam2=chebev(-1.0,1.0,c2, NUSE2,xx); 
*gampl= *gam2-x*(*gam1); 

*oammi= *gam2+x*(*gam1); 


i 
double chebev(double a, double b, double c[], int m, double x) 
{ 
void nrerror(char error_text[]); 
double d = 0.0, dd = 0.0, sv, y, y2; 
int J; 
if ((x-a)*(x-b) > 0.0) nrerror("x not in range in routine chebev"); 
y2 = 2.0*(y=(2.0*x-a-b)/(b-a)); 
for § = m - 1; j >= 1; j--) { 
sv =d; 
d = y2*d - dd + c[j]; 
dd = sv; 
} 
return y*d - dd + 0.5*c[0]; 
} 
#undef EPS 
#undef FPMIN 
#undef MAXIT 
#undef XMIN 
#undef PI 
#undef NRANSI 
#undef NUSE1 
#undef NUSE2 


#undef RTPIO2 
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APPENDIX J. ANSI C CODE FOR LEVIN ALGORITHM OPERATING ON 
EXTINCTION EFFICIENCY FACTOR 


/* Thesis Program for applying the Levin method on the extinction efficiency factor and determining */ 
/* the resulting rate of convergence */ 

/* LT Brian Johnson */ 

/* Compiler: Borland C++ Ver. 5.0 */ 

/* File Name: mielevin.c */ 


#include <stdio.h> 
#include <stdlib.h> 
#include <stddef.h> 
#include <math.h> 
#include “nrutil.h" 
#include “nrutil.c" 

#define NRANSI 

#define EPS 1.0e-10 
#define FPMIN 1.0e-30 
#define MAXIT 10000 
#define XMIN 2.0 

#define PI 3.141592653589793 
#define RTPIO2 1.2533141 
#define NUSE1 5 

#define NUSE2 5 

#define MAX 1036 


struct bessVec { 


double rj, ry, rjp, ryp; 
ie 


struct complex { 
double real, imag; 


}; 


void assignBessels(struct bess Vec* besselsPtr, double rj)dummy, double rydummy, double rjpdummy, 
double rypdummy) 
{ 


besselsPtr -> rj) = rjdummy; 
besselsPtr -> ry = rydummy; 
besselsPtr -> rjp = ryjpdummy; 
besselsPtr -> ryp = rypdummy; 


} 

struct bessVec sphbes(double, int); /* Function Prototype */ 

struct bessVec bessjy(double, double); /* Function Prototype */ 
double levin(); /* Function Prototype */ 
double bico(int, int); /* Function Prototype */ 
double factIn(int); /* Function Prototype */ 
double gammIn(double); /* Function Prototype */ 


struct bessVec sphBessels_n; 
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struct bessVec bessels_nPlusHalf; 

struct bessVec bessels_ nMinusHalf; 

struct complex a_n, b_n; /* Mie coefficients */ 
struct complex zeta_n; 

struct complex zeta_nMinus1; 


struct complex denom; /* denominator in Mie coefficients */ 
double x = 1000.0; /* size parameter */ 
double m = 2.0; /* refractive index (non-absorbing spheres) */ 


double J. nMinusHalf, J_ nPlusHalf; 
double A_ factor; 


double psi_nMinus1; /* Ricatti-Bessel fn 1st kind, initialized */ 
double psi_n; 

double chi_n; /* Ricatti-Bessel fn 2nd kind */ 
double nOverx; 

double aBraces, bBraces; 

double numer; /* numerator in Mie coefficients */ 
double realDenom; /* rationalized denominator */ 
double sum_QextOld, sum Qext, Q ext; /* Extinction efficiency factor */ 
double TIMAX+1], tt MAX+1]; 

double a; 

int n= 0; /* loop counter, also order of bessel function */ 

int mArray[] = {1.25, 1.33, 1.44, 1.50, 1.55, 2.0}: 

int k; 

main() 


{ 

T[O] = t[0] = 0.0; 

n= 0; 

printf("\n\ni T 1 t1 Pfi]/Qfij\n"); 


bessels_nPlusHalf = bessjy(m*x, 0.5); /* initialize for A_factor */ 
J_nPlusHalf = bessels_nPlusHalf.rj; 


psi_n = sin(x); 
zeta_n.real = psi_n; /* initializing zeta */ 
zeta_n.imag = cos(x); 


sum Qext= 0.0; 
while (n<=MAX) 
{ 

n=ntl; 


J_nMinusHalf = J_nPlusHalf; 
psi_nMinus1 = psi_n; 
zeta_nMinus1 = zeta_n; 
sum_QextOld = sum_Qext; 


bessels_nPlusHalf = bessjy(m*x, n + 0.5); /* calculating A factor */ 
J_nPlusHalf = bessels_nPlusHalf.1; 
A_factor = (J_nMinusHalf/J_nPlusHalf) - (n/(m*x)); 
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sphBessels_n = sphbes(x, n); /* j's and y's for Ricatti Bess fn's */ 
psi_n = x * sphBessels_n.xj; 

chi_n = -x * sphBessels_n.ry; 

zeta_n.real = psi_n; 

zeta_n.imag = chi_n; 


nOverx = n/x; 
aBraces = (A_factor / m + nOverx); 
bBraces = (A_factor * m + nOverx); 


numer = aBraces*psi_n - psi_nMinus1; 
denom.real = (aBraces * zeta_n.real) - zeta_nMinus1.real; 


denom.imag = (aBraces * zeta_n.imag) - zeta_nMinus]1.imag; 
realDenom = (denom.real*denom.real + denom.imag*denom.imag); 


a_n.real = (numer * denom.real)/realDenom; 
a_n.imag = (numer * denom.imag)/realDenom; 


numer = bBraces*psi_n - psi_nMinus]; 
denom.real = (bBraces * zeta_n.real) - zeta_nMinus1.real; 


denom.imag = (bBraces * zeta_n.imag) - zeta_nMinus].imag; 
realDenom = (denom.real*denom.real + denom.imag*denom.imag); 


b_n.real = (numer * denom.real)/realDenom; 
b_n.imag = (numer * denom.imag)/realDenom; 


t[n] = 2.0/(x*x)*(2*n + 1)*(a_n.real + b_n.real); 
T[n] = T[n-1] + t[n]; 


} /* end while loop */ 
a = levinQ; 
printf("\nx = %3.2f m = %1.2f \n\n", x, m); 


return 0; 
} 


struct bessVec sphbes(double x, int n) 

{ 
struct bess Vec tempsphbessels; 
struct bessVec bessels; 


double sj, sy, Sjp, Syp; 
void nrerror(char error_text[]}); 


double factor, order; 


if (n < 0 || x <= 0.0) nrerror("bad arguments in sphbes"); 
order=n+0.5; 


bessels = bessjy(x,order); 
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factor=RTPIO2/sqrt(x); 
sj=factor*bessels.rj; 
sy=factor*bessels. ry; 
sjp=factor*bessels.rjp-(sj)/(2.0*x); 
syp=factor*bessels. ryp-(sy)/(2.0*x); 


assignBessels(&tempsphbessels, sj, sy, sjp, Syp); 


return tempsphbessels; 


struct bess Vec bessjy(double x, double xnu) 


{ 


void nrerror(char error_text[]}); 


struct bess Vec tempBessels; 


double 1j, ry, tip, typ; 


void beschb(double x, double *gam1, double *gam2, double *gampl, 
double *gammi); 

int 1,1sign,1nl; 

double a,b, br,bi,c,cr,ci,d,del,del 1 ,den,di,dir,dli,dr,e,f,fact,fact2, 
fact3 ff, gam,gam],gam2,gammi, gampl,h,p,pimu,pimu2,q.r.1j1, 
rll ,xjmu,rp1 rjpl,rjtemp,ry 1 ,rymu,rymup,rytemp,sum,sum 1, 
temp, w,x2,xi,xi12,xmu,xmu2; 


if (x <= 0.0 || xnu < 0.0) nrerror("bad arguments in bessjy"); 
nl=(x < XMIN ? (int)(xnu+0.5) : IMAX(0,(int)(xnu-x+1.5))); 
xmu=xnu-nl; 
xmu2=xmu*xmu; 
xi=1.0/x; 
xi2=2.0* x1; 
w=xi2/PI; 
isign=1- 
h=xnu*¥ x1; 
if (h < FPMIN) h=FPMIN; 
b=xi2* xnu; 
d=0.0; 
c=h; 
for (i=1;1<=MAXIT;i++) { 
b += x12; 
d=b-d; 
if (fabs(d) < FPMIN) d=FPMIN; 
c=b-1.0/c; 
if (fabs(c) < FPMIN) c=FPMIN; 
d=1.0/d; 
del=c*d: 
h=del*h; 
if (d < 0.0) isign = -isign; 
if (fabs(del-1.0) < EPS) break; 
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if (1 > MAXIT) nrerror("x too large in bessjy; try asymptotic expansion"): 
rl=isign*FPMIN; 
npi=h* nl; 
nli=rl; 
npl=npl; 
fact=xnu*xi; 
for (=nl;1>=1;1--) { 
rjtemp=fact* ql+rypl; 
fact = x1; 
rjpl=fact* rjtemp-rl; 
rl=rjtemp; 
j 
if (rl == 0.0) ml=EPS; 
f=rjpl/rjl; 
if (x < XMIN) { 
x2=0.5*x; 
pimu=PI*xmu; 
fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimu)); 
d= -log(x2); 
e=xmu*d; 
fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e); 
beschb(xmu,&gam1,&gam2,&gampl,&gammi); 
ff=2.0/PI*fact*(gam 1*cosh(e)+gam2* fact2 *d); 
e=exp(e); 
p=e/(gampl*PI); 
q=1.0/(e*PI* gammi); 
pimu2=0.5*pimu; 
fact3 = (fabs(pimu2) < EPS ? 1.0 : sin(pimu2)/pimu2); 
r=PI* pimu2* fact3 *fact3; 
c=1.0; 
d = -x2* x2; 
sum=ff+r*q; 
sum 1=p; 
for (i=1;i<=MAXIT;i++) { 
ff=(1*ff+ p+q)/(i*i-xmu2); 
c *= (d/i); 
p = (-xmu); 
q /= (t+xmu); 
del=c*(ff+r*q); 
sum += del; 
del 1=c*p-i*del; 
sum! += dell; 
if (fabs(del) < (1.0+fabs(sum))*EPS) break; 
5 
if (i > MAXIT) nrerror("bessy series failed to converge"): 
rymu = -sum; 
ry! = -sum1*x12; 
rymup=xmu*xi*rymu-ry 1; 
rjmu=w/(rymup-f*rymu); 
} else { 
a—0.25-xmu2; 
p = -0.5*xi; 
q=1.0; 
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br=2.0* x; 
bi=2.0; 
fact=a*xi/(p*ptq*q); 
cr=br+q* fact; 
ci=bit p* fact; 
den=br*br+bi*b1; 
dr=br/den; 
di = -bi/den; 
dlr=cr*dr-ci* di; 
dli=cr*di+ci*dr; 
temp=p* dir-q*dli; 
q=p*dlit+q*dlr; 
p=temp; 
for (i=2;1<=MAXIT;i++) { 
a += 2*(1-1); 
bi += 2.0; 
dr=a* dr+br; 
di=a*ditbi; 
if (fabs(dr)+fabs(di) < FPMIN) dr=FPMIN; 
fact=a/(cr*cr+ci*ci); 
cr=br+cr* fact; 
ci=bi-ci* fact; 
if (fabs(cr)+fabs(c1) < FPMIN) cr=FPMIN; 
den=dr* dr+di*di; 
dr /= den; — 
di /= -den; 
dilr=cr* dr-ci* di; 
dli=cr* di+ci*dr; 
temp=p*dlr-q*dli; 
q=p*dlit+q*dlr; 
p=temp; 
if (fabs(dlr-1.0)+fabs(dli) < EPS) break; 
} 
if (i > MAXIT) nrerror("cf2 failed in bessjy"); 
gam=(p-f)/q; 
rmu=sqri(w/((p-f)*gam+q)); 
rymu=SIGN(jmu,yy)); 
rymu=rjmu* gam, 
rymup=rymu*(p+q/gam); 
ry 1=xmu* xi*rymu-rymup; 
} 
fact=rjmu/zj]; 
n=l *fact; 
p=np1* fact; 
for (i=1;i<=nl;i++) { 
rytemp=(xmuti)*xi2* ry 1-rymu; 


rymu=ry1; 
ry 1=rytemp; 
} 
ry—rymu; 


ryp=xnu*xi*rymu-ry 1; 


assignBessels(&tempBessels, rj, ry, Hp, Typ); 
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return tempBessels; 
j 


void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi) 
t 
double chebev(double a, double b, double c[{], int m, double x); 
double xx; 
Static double cl[] = { 
-1.142022680371172e0,6.516511267076e-3, 
3.087090 17308e-4,-3.470626964e-6,6.943764e-9, 
3.6780e-11,-1.36e-13}; 
static double c2[] = £ 
1.843740587300906e0,-0.076852840844786e0, 
1.271927136655e-3,-4.971736704e-6,-3.3 126120e-8, 
2.42310e-10,-1.70e-13,-1.0e-15}; 


xx=8.0*x*x-1.0; 

* 9am 1=chebev(-1.0,1.0,c1, NUSE1,xx); 
* 9am2=chebev(-1.0,1.0,c2, NUSE2,xx); 
*gampl= *gam2-x*(*gam1); 

*gammi= *gam2+x*(*gam1); 


} 
double chebev(double a, double b, double c[], int m, double x) 
{ 
void nrerror(char error_text[]); 
double d = 0.0, dd = 0.0, sv, y, y2; 
int J; 
if ((x-a)*(x-b) > 0.0) nrerror("x not in range in routine chebev"); 
y2 = 2.0*(y=(2.0*x-a-b)/(b-a)); 
lent m=); >= 1; jm 
sv=d; 
d = y2*d - dd + cf]; 
dd = sv; 
return y*d - dd + 0.5*c[0]; 
i 
double levin() 


{ 
double PIMAX+1], QIMAX+1]; 
double commonTerm; 
int N; 
P[O}] = Q(0] = 0.0; 


for (N = 1; N <= MAX; N++) 


for (k= 1; k <= N; k++) 
{ 
commonTerm = pow(-1,k)*pow(k,N-1)/t[k] *bico(N,k); 
P[k] = P[k-1] + commonTerm*T[k]; 
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Qik] = Q[k-1] + commonTerm; 
} 
printf("\n%d %f Nw %E",k- 1, TIN]. t{N], PIN]/Q{N]); 


} 

printf("\n\n"); 

return (P[MAX]/Q[MAX]); 
} 


double bico(int n, int k) 
{ double factln(int n); 
return floor(0.5+exp(factln(n)-factin(k)-factln(n-k))); 


; 
double factIn(int n) 
{ double gammiIn(double xx); 
static double a[101]; 
if (n <= 1) return 0.0; 
if (n <= 100) return a[n] ? a[n] : (a[n]=gammln(n+1.0)); 
else return gammln(n+1.0); 
j 
double gammin(double xx) 
{ double x,y,tmp,ser; 
static double cof[6]={76. 18009 172947146,-86.5053203294 1677, 
24.01409824083091,-1.23 1739572450155, 
0.1208650973866179e-2,-0.5395239384953e-5}- 
int j; 
Y—X= XX; 
tmp=x+5.5; 
tmp = (x+0.5)*log(tmp); 
ser=1.000000000190015; 
for §=0:j<=5;j++) ser += coffj]/+y; 
return -tmptlog(2.5066282746310005*ser/x); 
} 
#undef EPS 
#undef FPMIN 
#undef MAXIT 
#undef XMIN 
#undef PI 
#undef NRANSI 
#undef NUSE1 
#undef NUSE2 
#undef RTPIO2 
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